####复杂热图###########################
rm(list = ls())
options(stringsAsFactors = F)

library(clusterProfiler)
library(msigdbr)
library(GSVA)
library(GSEABase)
library(pheatmap)
library(limma)
library(colorRamps)
library(ComplexHeatmap)
library(vegan)
library(circlize)
library(Cairo)
library(corrplot)

setwd('父级路径')
##提取自己的样本
FPKM=read.table("convert_exp.txt",sep="\t",header=T,check.names=F, row.names = 1)

stadns<-'ista'

if (stadns=='ista') {
  message("标准化了FPKM")
  dat <- as.matrix(log2(FPKM+1))
}else{
  dat <- as.matrix(FPKM)
}

##这一步是提取自己的基因集
genelist <- read.table("工作簿3.csv",sep=",",header=T,check.names=F)
geneset <- lapply(genelist, function(col) col[!is.na(col) & col != ""])
huangPar <- gsvaParam(exprData = dat, geneSets = geneset,kcdf = "Gaussian",minSize  =1,maxSize = )
gsva_mat <- gsva(huangPar)#调用所有核
counts=read.table("GSE126044_counts.txt",sep="\t",header=T,check.names=F, row.names = 1)

if (stadns=='ista') {
  message("标准化了counts")
  dat <- as.matrix(log2(counts+1))
}else{
  dat <- as.matrix(counts)
}


genelist <- read.table("工作簿3.csv",sep=",",header=T,check.names=F)
geneset <- lapply(genelist, function(col) col[!is.na(col) & col != ""])

huangPar2 <- gsvaParam(exprData = dat, geneSets = geneset,kcdf = "Gaussian")
gsva_mat2 <- gsva(huangPar2)#调用所有核
merge <- cbind(gsva_mat2,gsva_mat)
write.csv(gsva_mat,"gsva_go_matrix.csv")

# 生成图
cluser = read.csv(file = 'ConsensusClusterResults.k=3.consensusClass.csv', sep = ',', header = F)
rt <- read.csv(file = 'gsva_go_matrix.csv', sep = ',', header = T, row.names = 1)
rownames(cluser) <- cluser[,1]
names(cluser)[2] <- "Cluster"
cluster1_rows <- rownames(cluser[cluser$Cluster == "1", ])#看一下1有多少，这里是221
rt_intersect1 <- rt[, intersect(colnames(rt), cluster1_rows)]
##再来一遍
cluster2_rows <- rownames(cluser[cluser$Cluster == "2", ])#看一下2有多少，这里是276
rt_intersect2 <- rt[, intersect(colnames(rt), cluster2_rows)]
cluster3_rows <- rownames(cluser[cluser$Cluster == "3", ])#看一下3有多少，这里是176
rt_intersect3 <- rt[, intersect(colnames(rt), cluster3_rows)]

merged <- cbind(rt_intersect1, rt_intersect2, rt_intersect3)
merged<-as.matrix(merged)
merged <- decostand(merged,"standardize",MARGIN = 1)#标准化
col_fun = colorRamp2(c(-2, 0, 2), c("#1f60a4","white","#a40d20"))#,space = "RGB"
col_Cluster1 <- c("1" = "#B49FDA", "2" = "#72AFD3", "3" = "#5FDFFF")
cluser <- cluser[colnames(merged),]
column_ann <- HeatmapAnnotation(Cluster=cluser$Cluster,
                                col = list(Cluster = col_Cluster1)
)

png(file="p1.png",width=980, height=980)

pictureTemp

#  Heatmap(merged,
#               col = col_fun,
#               cluster_rows = T,##打开行聚类
#               cluster_columns = T,##关闭列聚类##show_row_dend = T, show_column_dend = F聚类但是不显示聚类树
#               show_row_names = T,##行名
#               show_column_names = T, ##列名
#               #column_names_rot = 2,##列名旋转方向
#               column_names_side = "top",#上方列名
#               #column_dend_side = "bottom"，#同上
#               ##row_km = 2, #k-mean法分割行
#               ##column_split = split,##自定义分割
#               top_annotation = column_ann
# )
dev.off()
